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We present a brief overview of the methods for making statistical inference (testing statistical hypotheses, construction of confi- 
dence and/or prediction intervals and regions) about linear functions of the fixed effects and/or about the fixed and random effects 
simultaneously, in conventional simple linear mixed model. The presented approach is based on solutions from the Henderson's 
mixed model equations. 
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The applications of data analysis based on the statistical lin- 
ear mixed model, as a natural generalization of the a naly sis of 
variance methods and the ANOVA models, (see e.g. 14411 . H 1 511 . 
l36ll ). are widespread. Such applications with analytical meth- 
ods based on linear mixed models include different fields of the 



biomedical and technical research, (see 115 611 and/or 111 111 ). For 
illustration, here we shall mention just few of them: e.g . ge- 
netics with its microarray experiments, |0K llH, 10, |74|, the 
plant and animal breeding in agricultural, [5], statistical meta- 
analysis in medical research, 111 811 . neurophysiology, 115 ill , as 
well as different technical applications, like e.g. calibration of 
devices, derivation of the tolerance intervals for industrial appli- 
cations, interlaboratory comparisons in metrology, and methods 
for expression the uncertainties in measurements, see e.g.Jfj], 

El, S, Jp, JH, 0, m, HI, H, m, & m, 
izaundHT 

Although the linear mixed models and the methods for sta- 
tistical inference based on such models have been recognized 
and used for long time by the researchers in different fields, 
it seems that some sort of misunderstanding of the principles 
and/or the technical details (of the used methods for statisti- 
cal inference based on such linear mixed models) may lead to 
improper usage of the implemented methods and algorithms. 
Moreover, there are still some further open theoretical prob- 
lems (like e.g. methods for testing and constructing confidence 
intervals/regions about the variance components, see e.g. 0, 

H, 0i, HI, 0, 0, & & S, mil). 

So, the main goal of the paper is to present a brief overview 
of the standard (conventionally used) methods for making sta- 



tistical inference (in particular the methods for testing statistical 
hypotheses and the methods for construction of the confidence 
and/or prediction intervals/regions) about linear functions of the 
fixed effects and/or about the fixed and random effects simul- 
taneously, in conventional simple linear mixed model, (with 
pointing to potential problems which may appear based on us- 
age of these methods), and to present some of the recently de- 
veloped improvements, as well as some generalizations, together 
with relatively detailed technical description of the model and 
the methods. The presented approach is based on the elements 
of the solution of the Henderson's mixed model equations. 

2. Henderson's mixed model equations 

We consider the linear mixed model (LMM) in the follow- 
ing form 



y = Xb + Zu + . 



(1) 
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with y being an 71-dimensional vector of observations, b be- 
ing the p-vector of fixed effects, u being the r-vector of ran- 
dom effects with E(u) = and Var(u) = G, and e being the 
n-vector of random (measurement) errors with E(e) = and 
Var(e) = R, where R is assumed to be strictly positive definite 
variance-covariance matrix of e. The (« x /;>)-matrix X and the 
(n x r)-matrix Z are the known design matrices. Typically, we 
can write Zu = Yn=\ Zi u i> where the (n x r,) matrices Z, and the 
/•/-dimensional random effects m,, i = l,.,.,s, could be speci- 
fied from the structure of the model. 

The main goal of this paper is to present an overview of the 
methods for making statistical inference about linear functions 
of the fixed effects b and the random effects u, i.e. about K'b 
and/or about w = A'(b',u')' = K'b + L'u for given (suitable) 
coefficient matrices A, resp. K and L. 
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Henderson in 12311 developed a set of equations, termed as 
the mixed model equations (MMEs), that simultaneously yield 
the best linear unbiased estimator (BLUE) of Xb (or any vec- 
tor of estimable linear functions K'b) and the best linear un- 
biased predictor (BLUP) of u (or any vector w = K'b + L'u, 
provided K'b is estimable), under the assumption that the co- 
variance structure is known. 

The MMEs were derived based on the normality assump- 
tions, i.e. u ~ N(Q,G), e ~ N(0,R), with Cov(u,e) = 0, for 
known variance-covariance matrices G and R. Thus, the joint 
probability density function (pdf) of the random vector (y', u')' 
is given as 



f(y,u) = f(y\u)f(u) 

= — — exp 1 — (y - Xb - Zu)'R~\y - Xb - Zu) 



1 1,, 

! — exp< — u G u 

(2jtYI 2 \G\ 1 I 2 F \ 2 



By solving the ML equations for b and u, i.e. 

dfiy, u) = df(y, u) _ 

db du 
we get the MMEs in the following form 



X'R X X X'R X Z 
Z'R l X Z'R l Z + G l 



X'R~ l y 
Z'R l y 



(2) 



(3) 



(4) 



The left-hand side matrix of © will be termed as the Hender- 
son's MME matrix, here denoted by H, i.e. 



H = (X,Z)'R-\X,Z) + (OJrYG-'iOJr), 



>n-h 



(5) 



where by we denote a zero matrix with suitable dimensions, 
here (r x p). Alternatively, 



X'R l X X'R X ZG 
Z'R l X W~ l 



X'R l y 
Z'R l y 



(6) 



where W = (/ + Z'R~ l ZG)~ l . Notice, that based on ©, there 
is no need to restrict the variance-covariance matrix G to be 
strictly positive definite. This version of MMEs is preferred for 
numerical evaluations, if G can be a bad conditioned matrix. 

Given the variance-covariance matrices G and R, let us de- 
note as C the following matrix of coefficients 



C = 



Cn C12 

C<1\ C22 

X'R l X X'R l Z 
Z'R l X Z'R l Z + G l 
I n \/ X'R l X X'R l ZG 



G J\ Z'R l X 



(7) 



where by A~ we denote any g-inverse of the matrix A. 

Let b and 5 be any solution to the MMEs |@}. Notice that 
based on b and v, the solutions from ©, we can reconstruct 5 
by u — Gv. Then the BLUE of the vector of linear estimable 



functions of the fixed effects K'b, see e.g. 04911 . is 



BLUE(K'b) = K' (X'V~ 1 X\ X'V~ l y = K'b, 



(8) 



where K' is a (q x p)-matrix of coefficients of the estimable 
linear function K'b, i.e. K = X'A for some matrix A, and V = 
Z'GZ + R. The BLUP of the vector of linear functions of the 
fixed and random effects, say K'b + L'u, is 



BLUP(K'b + L'u) 



= BLUE(K'b) 

+L'GZ'V~ l (y - BLUE(Xb)), 
= K'b + L'u, (9) 



where L' is an arbitrary (gxr)-matrix of coefficients, and BLUE(Xb) 
Xb. 

Important properties of the solutions of the MMEs are sum- 
marized bellow, for more details see e.g. 113 811 : 



2. 
3. 
4. 
5. 



1 . In the class of linear unbiased predictors, BLUP maxi- 
mizes the correlation between u and S. 
K'b is BLUE of the set of estimable linear functions K'b. 
E (u\u) = u. 
11 is unique. 

K'b + L'u. is BLUP of K'b + L'u provided that K'b is 
estimable. 

6. VariK'b) = K'C U K. 

7. VariK'b + L'u) = K'C U K + L'(G - C 2 i)L. 

8. VariK'b + L'u) - {K'b + L'u)) = (K' , L')C(K' , L')' . 

9. Cov(K'~b,u'\ = 0. 

10. Cov(K'b,u') = -K'C n - 

11. Cov (K'b,u' -u') = -K'Cn- 

12. Var(u) — Cov(u, u') — G — Cn- 

13. Var (u — u) = Cn- 

In this paper we shall consider only a special form of the 
model (Q]l — a conventional simple LMM with normally dis- 
tributed errors and random effects. That is, we shall assume mu- 
tually uncorrected (independent) normally distributed random 
effects u 1 , . . . , u s and e with E(uf) = for i = 1 , . . . , r, E{e) = 0, 
Cov(uj, uj) = for i + j, and Cov(u t , e) = for all i = 1, . . ., s. 
Further, we shall assume Var(ui) = ajl n , i = l,...,s, with 
r = r /> an d Var(e) = cr 2 s+l I„. Hence, 



E(y) = Xb, and Var(y) = j] ojZfi + cr 2 s+l I„, 



(10) 



with cr 2 = (cr 2 , . . ., <x 2 , <x 2 +1 ) being the vector of variance com- 
ponents with the parameter space specified by cr 2 > for i = 
1, . . . , s, and cr 2 +l > 0. However, in order to avoid possible 
technical and numerical problems, it is reasonable to assume 
that the true parameter cr 2 = (cr 2 , . . . , cr 2 , cr 2 +l ) is in the in- 
terior of this parameter space. So, here we shall assume that 
cr 2 > for i = 1, . . . , s + 1, 

In other words, we shall assume y ~ N(Xb, V), with V = 
Var(y) = ZGZ' + R, where G is (r x r) diagonal matrix, G 
Vaii u ) = diag(cr 2 / r ), and R is (n x n) diagonal matrix, R = 
Var(e) = cr 2 +1 I„, with cr 2 > for i = 1, . . . , s + 1. 

If the variance components cr 2 = (cr 2 , . . . , cr 2 , cr 2 +l ) are un- 
known, they can be (and in general must be) estimated from 
the observed data by any reasonably effective and computation- 
ally efficient method, like e.g. by the methods based on mo- 
ments (the minimum variance (norm) quadratic estimation) or 
the methods based on likelihood function (ML or REML). 
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There are several efficient implementations for estimation 
of the variance components in general LMMs. One method 
used to fit such LMMs is the expectation-maximization (EM) 
algorithm, see [34], where the variance components are treated 
as unobserved nuisance parameters in the joint likelihood. Cur- 
rently, such methods are implemented in the major statistical 
software packages SAS (Proc MIXED) and R (lme in the nlme 
library). In particular, Proc MIXED uses a ridge-stabilized Newton- 
Raphson algorithm to optimize either a full (ML) or residual 
(REML) likelihood function, see also iS. j3al. lieoTl. and lf40Tl. 

However, here we present a relatively simple method, based 
on repeated iterative solving of the MMEs, suggested by Searle, 
Casella and McCulloch in fl49[] . The elements of MMEs are 
used for setting up iterative procedures for simultaneous es- 
timation of the variance components cr 2 , . . . , <j 2 s , <x 2 +l and the 
empirical versions of the BLUE of b and the BLUP of u, in the 
simple LMM ([T0>. 

The algorithm provides solution to the maximum likelihood 
(ML) or the restricted maximum likelihood (REML) equations 
for estimating variance components, see e.g. |17], |39], lfl9ll . 
1 321. and l49ll . The algorithm canbealsoused for estimation of 
the related Fisher information matrices for ML and/or REML 
estimators of the variance components (i.e. the inverse of the 
asymptotic variance-covariance matrix of the ML/REML esti- 
mators). Moreover, it can be also used for computing the min- 
imum norm quadratic estimates MINQE(I) (realizations of the 
invariant minimum norm quadratic estimators) or the MINQE(U,I) 
(invariant and unbiased minimum norm quadratic estimators) 
of the variance components, for more details see e.g. 1I33I1 . |H2l, 
and H. 

The final solutions of such iterative procedure will be de- 
noted by b, u = (u[, . . ., u' s )', and & 2 = (& 2 , . . . , <x 2 +I )'. Sim- 
ilarly, we shall use the adequate notation G, R, and C for the 
estimated versions of matrices G, R, and C. The solutions b 
and u satisfy the MMEs (O if the unknown matrices G and 
R are replaced by the estimated versions G and R. Finally, 
based on cr 2 , the important output of the algorithm is the esti- 
mated Fisher information matrix, say /ml(<5" 2 ) or Ireml(ct 2 ), re- 
spectively. Consequently, it provides the estimated asymptotic 
variance-covariance matrix of the estimated variance compo- 
nents cr 2 , say £ = (/mz.(<5~ 2 )) or £ = (i RE ml{^ 2 )) , provided 
that the inverses do exist. For detailed description of the algo- 
rithm see Section [Appendix ~B~| 



3. Standard methods for statistical inference on fixed and 
random effects 

Here we consider the problem of making statistical infer- 
ence about q linear functions of the fixed effects b and the ran- 
dom effects u, i.e. about A' (b', u'Y - K'b + L'u where A is 
((p + r) x g)-dimensional full-ranked matrix with estimable K'b 
(i.e. K - X'A for some matrix A). 

Let b and u are the solutions of the MMEs (O, so w = 
A' (b', u')' - K'b + L'u is the best linear unbiased predictor 
(BLUP) of w = K'b + L'u. Then, according to the properties 6 



and 8 of Section [5] the variance of K'b and the mean squared 
error (MSE) of w are given by 



and 



Var(K'b) = K'C U K, 



MSE (w) = E((w- w) (w - w)') 

= Variw- w) = A'CA = Af*. 



(11) 



(12) 



Notice that the MSE matrix of w, M vT ,, functionally depends on 
the variance components cr 2 = (cr 2 , . . . , <x 2 , <x 2 +1 ) . 

If the variance components <x 2 = (cr 2 , . . . , cr 2 , <x 2 +1 ) are 
known, based on the model assumptions and from dTTb and 
(fT2l . we trivially get the pivot, Wald-type statistic, useful for 
making statistical inference about K'b (e.g. testing a null hy- 
pothesis Hq : K'b = K'bo for some bo) and/or about the variable 
w = K'b + L'u with their exact (null) distribution: 



Q = (K'b - K'bo)' (K'CuK)- 1 (K'b - K'bo) ~ x], (13) 



and 



Q = (w-wY (A'CA)- 1 (w-w)~xi 



(14) 



where x\ denotes the chi-squared distribution with q = rank(A"') = 
rank(A') degrees of freedom. 

If the variance components are unknown and the estimated 
values <r 2 = (cr 2 , . . . ,fr 2 +l ) are available together with C, a 
commonly used test statistic for fixed effects hypothesis Hq : 
K'b = K'b , is based on K'b and C n : 

F = - (K'b - K'bo)' (K'C n K) 1 (K'b - K'bo) , (15) 

where K'b denotes the empirical version of the best linear un- 
biased estimator K'b of K'b (i.e. version with the estimated 
variance-covariance components). Notice that Cn =\X'V~ 1 X) , 

see e.g. [49] (Eqn. (55) p. 276), and consequently C n = (x'V~ l x)~ 
where V = ZGZ' + R. 

As a generalization, for making simultaneous statistical in- 
ference on the fixed as well as the random effects, i.e. on w — 
A' (b' , u')' (e.g. construction of the prediction region) based on 
the empirical BLUP (EBLUP), i.e. the predictor w = A' (b', u')' 

(where b and u are solutions of the MMEs with estimated R and 
G), it is natural to consider the following statistic 



1 / - \-i 

F = - (vv - w) I A CA) (w-w). 
q V ) 



(16) 



where q is rank of the matrix A'. 

As a special case, if w is a one-dimensional function given 
by w = X (b', u')' - k'b + l'u, in analogy with ( fTBT ) and ( fTBT ), it 
is natural to consider the pivot statistic 



t - 



k'b - k'bo 



y/k'Cuk 
and/or its generalization 
w-w 
y/xCA 



(17) 



(18) 



where w = A' (h',u'Y is the EBLUP of w. 

The (null) distribution of the statistics ([T71 > and (TT~8T > is com- 
monly approximated by the Student's f-distribution with v de- 
grees of freedom (DF), estimated by applying the Satterthwaite's 
approximation. The (null) distribution of the statistics (TT~5T > and 
( TToT ) is commonly approximated by the Fisher-Snedecor's F- 
distribution with v\ and V2 degrees of freedom, where v\ — q 
and V2, the denominator degrees of freedom (DDF), where V2 
is typically estimated by a generalization of the Satterthwaite's 
method, as suggested e.g. by Fai and Cornelius in 11311 . or al- 
ternatively, by applying moment based approximation for the 
F-distribution. The explicit expressions for DF and DDF esti- 
mators of ( fTTl i. ( TT~8T >, (TT3T > and (TToT l are given in Sections [XJJand 
331 



3.1. DF estimated by the Satterthwaite 's method 



Giesbrecht and Bums in [ 16], (see also [|37|], llj], and [50]), 
suggested to approximate the null distribution of the pivotal 
quantity dTvT > by the Student's f-distribution with v degrees of 
freedom (DF), where v is the Satterthwaite's approximatiorQ of 
the (unknown) v, see 114611 . IHTll . i.e. 



t = 



k'b - k'b 



with 



v k = 



y/k'Cnk 
2(k'C n kf 



2{k'C u kf 



Var (k'Cuk) gfigk 



(19) 



(20) 



where Var(k'C\\k) denotes the estimated value of Var (k'C\\k). 

The suggested estimator of Var(k'C\\k) = g' k f.gk is based 
on the estimated version of the Taylor series expansion of the 
variance of the estimator k'b (BLUE), i.e. Var [k'b) = k'C u k, 
with respect to the variance components <x 2 = (<x 2 , . . . , <x 2 , cr 2 +1 ). 
Here, £ is the estimated (asymptotic) variance-covariance ma- 
trix of the estimators (e.g. REML estimators) of the variance 
components a 2 , and gk is the estimated version (evaluated at 
the estimated values of the variance components <x 2 ) of the gra- 
dient gk of k'Cuk, with respect to the variance components a 2 , 
i.e. 



gk 



l d(k'C n k) 
do? 



S(.k'C n k) 

da] 
m'Cnk) 



(21) 



The Satterthwaite's approximation of the distribution of k'Cuk is based 
on assumption that v(t'Cni) /<x 2 ~ x~ for some parameters <x 2 and v. By 
comparing the first and the second moments of both random variables we get 
E{v {k'C\ \k) /a 2 ) = v and Var(y(k' Cnk) /cr 2 ) = 2v. From that we directly 

get a- 2 = E(k'C n k) and v = 2 (e (k'C n k)f I Var (k'C n k). As E(k' C n k) and 
Var(k'Cuk) depend on unknown parameters they should be estimated. So, we 
get the natural estimator as v = 2{k'C\\k^j IVar(k'C\ikj. 



As a generalization of the approach by Giesbrecht and Burns, 
it is natural to consider similar approximation for the distribu- 
tion of the pivotal quantity (flSt . i.e. 



(22) 



with 



va = 



y/A'CA 
2{A'CA) 2 2[A'CA) 2 



Var(A'CA) g' A I,gA 



(23) 



where g ^ is the estimated version of the gradient g^ of MSE (w) 
A'CA with respect to the variance components cr 2 , defined by 



gA = 




(24) 



For more details on computing gradients of the MSE(w) see 



Section Appendix A 



Provided that the estimated matrix C is available, e.g. as an 
output of the algorithm for estimating the variance components, 
the estimators gk and g,\ of the gradients (fJTJ and ( l24l i could 
be evaluated, by using the elements of the estimated matrix C 
(instead of C). 

For that, let us define A = CA and let A be decomposed into 
its subvectors such that A = (%'„,![,. ..,A.' S )', where Ao is p- 
dimensional subvector, and A/, i = 1, . . . , s, are r, -dimensional 
subvectors of A. Then, by using ( IA.3U from Section Appendix A. 3 
we get 



gA = 



■AM 



tA'HqA 



where Ho is given by 

H = {X,Z)'{X,Z) 



(25) 



X'X 
Z'X 



X'Z 
Z'Z 



(26) 



Consequently, as k'b is a special case of A' (b\ u')' = k'b + 
I'u with A = A(k) = (k' , Q' r )' , so we can use d25l > also for evalua- 
tion of gk by replacing A with \q = CA(k) ■ 

3.2. DDF estimated by the Fai-Cornelius method 

Fai and Cornelius in IU3I1 proposed a generalization of the 
Satterthwaite's method for multivariate linear functions of the 
fixed and random effects to approximate the (null) distribution 
of the statistic ( TTBT l by the Fisher-Snedecor F-distribution with 
v\ — q and V2 = v, i.e. with the estimated denominator degrees 
of freedom (DDF). 

As a straightforward generalization of the Fai-Cornelius ap- 
proach, it is natural to approximate the distribution of the F- 
statistic ( TTBT l. based on the multivariate function w = A' (b', u')' - 
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K'b + L'u and its empirical predictor vv = K'b + L'u, by the 
Fisher-Snedecor F-distribution with v\ — q and V2 — v degrees 
of freedom, where where 



v = 



2E 
E - q 



with 



g ~ X ^2 1{ * >2 > ' 



(27) 



(28) 



Here, lj.j denotes the indicator function and v,, for i = 1, . . . 
are the degrees of freedom, estimated by the Satterthwaite's 
method ( 1231 ). of the f-statistics ( fT8l for w, = A' t (£>', fi'j , where 
1;, i = 1, . . . , q, are the columns of the matrix Ape given by 



A 



FC 



AU, 



(29) 



and U denotes the unitary matrix of a spectral decomposition 
of a matrix A'CA, i.e. such matrix that U'A'CAU = S , where 
S is a diagonal matrix. 



, j+l j+l ~2 ~ 
'=1 7=1 • 7 



Then taking expectation of the square of the first-order term, 
and using the results in [28] and 12111 . we get the first-order 
approximation Mg$ of Msa, as 



M, 



, . dw ^dw' 
do~ 2 ' dcr 2 



E2>* 

i=l 7=1 
J+l J+l 

XX 2 '; C0V 

'=1 7=1 



' <9vv dw' 
deader 2 ) 

8 (w — w) d(w —w) 



dcr 2 



do- 2 



(33) 



where are elements of the variance-covariance matrix X of 
the estimator & 2 . 



For derivation of the approximation of Msw see Section Appendix A.' 
The second component of the EBLUP's MSE matrix Mga, in the 
simple LMM ( fTOb can be approximated by 



4. Statistical inference on fixed and random effects based 
on adjusted estimator of the MSE matrix of the EBLUP 



As argued by Harville in H22I1 . usage of the MSE matrix 
of the BLUP w, say Ma,, (or its estimated version, say Ma,), 
instead of the correct MSE matrix of the EBLUP vv, say Ma,, (or 
its estimated version, say M#), is inadequate, as the estimator 
Ma, = A'CA can severely underestimate the true MSE of the 
EBLUP w. As will be explained bellow, there are two main 
sources of such bias. For a comprehensive discussion on the 
problem and proposed solutions see also lE 7|. lE8ll . lEoll . 1251. 

liH, mi, M, m, JH, J3, d, Him and mr 



4.1. Decomposition of the EBLUP prediction error and its MSE 
The first source of the bias can be observed if we decompose 
the prediction error of the EBLUP vv. In particular, 



(vv - w) — (vv - w) + (vv - w) , 



(30) 



and consequently, based on unbiasedness of EBLUP and its in- 
dependence on BLUP, see |27], |28], ll20ll . and |21], we get the 
MSE matrix of vv in the form 



M & = My, + M 5& , 



(31) 



where Mg$ = E ((vv - vv) (vv - vv)') = Var (vv - vv), and thus, 
Ma, > Ma,. 

The MSE of the first component of the prediction error, Ma,, 
is given by (TT2l) . The MSE of the second component of the pre- 
diction error, Msa,, is not expressible in closed form, except for 
very simple special cases. Kackar and Harville in [28], see also 
i29fl and i30ll . suggested approximation of Mg$ based on first- 
order Taylor series approximation. In particular, a Taylor series 
expansion for vv - vv in o -2 = (cr 2 , . . . , fr 2 , 0" 2 +1 ) , as e.g. REML, 
about cr 2 = (cr 2 , . . ., cr 2 , cr 2 + ^j , gives approximation 

J+ 1 

(vv - vv) * (vv - vv) + 2^ ^3 ~ °f) 



M„ 



-XX--- 



(34) 



i=l 7=1 



where Cy, i, j = 1, ...,s+ 1, are given by ( 1A.40I I. or alterna- 
tively by 



j j+i j+i 

^* = -?XX^ M * J 



(35) 



i=l 7=1 



where the matrices M§ are given by ( 1A.321 I. ( IA.33b . ( IA.341 I. 
and dAJBl 

Consequently, we get the approximation M t - t . of the EBLUP's 
MSE matrix M,- v in the form 



M A = Mjj) + M s& 

j+i j+i 



J+l J+l 



(36) 



'=1 7=1 



where Ey are elements of the variance-covariance matrix of 
the REML estimator a 2 , and Mv 1 represent the second par- 
tial derivatives of the BLUP's MSE matrix Ma, with respect to 
the variance components cr 2 and cr 2 , i, j = 1, . . . , s+l, in simple 
LMM OB. 

4.2. Bias-corrected estimator of the EBLUP 's MSE matrix M v \, 

As the EBLUP's MSE matrix Ma,, as well as its approxi- 
mation M( v (which is a function of Z), depend on the unknown 
variance components cr 2 = (cr 2 , . . . , cr 2 ^ , for further appli- 
cations it is necessary to use its estimator, say M&. A natural 
option for such estimator would be 



M,-v = Ma, + M, 



(37) 



i.e. by using ( f36l >. where the true (unknown) vector of variance 
components a 2 is replaced by its estimator a 2 . Notice that £, 
the true variance-covariance matrix of the REML estimator a 2 
also depends on a 2 . So, the estimator (l37T > functionally depends 

on the elements of estimated variance-covariance matrix £. 

j— i 

Based on similar arguments as given by Alnosaier in yj] for 
the special case of empirical BLUE of the fixed effects, we can 

assume that M$k is approximately unbiased estimator of Nig®, 
for another formal justification see also 114 ill and 111 OB. 

However, as pointed out by Harville and Jeske in 112 ill . Prasad 
and Rao in [41], and in special case of fixed effects estimator by 
Kenward and Roger in M29I1 and 13011 . additional bias will appear 
if the estimator M vT , is used as an estimators of the MSE matrix 
Mf V in ( T37T i. In order to show that, let us expand M% in a 2 about 
a 2 , and then take expectation of this approximation, so 



E(M^)^M^ + Y j E(d 2 i -o 2 ) 



i=\ 
s+1 j+1 



Ite 2 



4zf>(H-^-"i))^ 



'=1 7=1 

.s+1 .s+1 



<0~ 



M« - M 6& , 



(38) 



where we have assumed that the first-order term could be ig- 
nored, and Mgcv is given by P5t . This could be informally jus- 
tified by the assumption that frj is approximately an unbiased 
estimator of a 2 , as was suggested in Il29fl . However, formal jus- 
tification was provided by Alnosaier in [ 1 ] and by Kenward and 
Roger in [30]. Kenward and Roger derived Taylor series ap- 
proximation for the bias of REML estimator, i.e. EUJ 2 - a 2 ^, 
and proved that in linear mixed models with linear parametriza- 
tion of the variance-covariance matrix V = Z'GZ+R, like e.g. in 
simple LMM dTOb . its first-order approximation is equal to zero. 

Hence, by combining d37b and d38b . we get the adjusted, 
bias-corrected estimator of the EBLUP's MSE matrix M*, given 
by 



M*_a = M* + 2Ma 



(39) 



The explicit form of the estimator (l39t in simple LMM (TTOb is 
given by ( IA.45b in Section Appendix A. 5 



4.3. Generalization of the Kenward-Roger method for statis- 
tical inference on fixed and random effects based on ad- 
justed estimator of the MSE matrix of the EBLUP 

For statistical inference about the vector of linear functions 
of fixed effects K'b based on its empirical BLUE, Kenward and 
Roger suggested in [29] to use the Wald-type statistic as a pivot, 
with adjusted covariance matrix of the empirical BLUE of the 
function K'b. 

Here we suggest to consider a generalization of the Kenward- 
Roger method for the inference about the vector of functions of 
fixed and random effects w = A'(b',u')' (which is useful for 
testing hypotheses about the fixed effects and for constructing 



the prediction regions for functions of the fixed and the random 
effects simultaneously), based on its EBLUP and the adjusted 
MSE matrix. For that we shall consider the Wald-type pivot 
F-statistic 



F = - (w - w)' {m^ (w-w), 



(40) 



where M^a is given by d39b . or (in its explicit form) by (IA.45b 



from Section Appendix A. 5 respectively. 

In accordance with [29] and [ 1], we suggest to approximate 
the (null) distribution of the scaled Wald-type F-statistic d40b 
by the Fisher-Snedecor ^-distribution with q and v degrees of 
freedom. In particular, 



approx. p 



q,v> 



(41) 



where the unknown parameters k and v should be estimated 
from the data. 

In analogy with derivation of the estimators presented by 
Alnosaier in flJJ] for the fixed effects problem, here we suggest 
the following estimators of the scale k and the denominator de- 
grees of freedom v: 



= 4 + 



E(v-2) 
2 + q 



where 



Q 
E 
V 
B 



qg-1 



V 
2E 2 ' 

i 

2 -(l + B), 
q v I 

2^(A, + 6A 2 ), 



(42) 



(43) 



and 



''=i ;'=i 

.S+1 J+1 

iz = ZZ £ y tr (^* li ^ )i ^ 1 ^)- (44) 



'=1 /=! 



By tr(A) we denote the trace of a matrix A, i.e. tr(A) = 2; Yjj 
M# = A'CA denotes the estimated version of My,, and Mi , 
i = 1, . . . , s+ 1, denote the estimated versions of the first partial 
derivatives of M s , defined by (IA.3 11 >. For more details and ex- 



plicit forms of the estimators A i and A2 see Section Appendix A. 6 
dA60b and iKM\ . 

In order to match the exact values for the scale k and the de- 
nominator degrees of freedom v for testing hypothesis on fixed 
effects in two special cases, in particular in the balanced one- 
way ANOVA and the Hotelling T 2 models, Kenward and Roger 
in 129ll suggested the modified estimators k* and v*, which can 



6 



be analogically generalized and used to approximate the (null) 
distribution of the scaled Wald-type F-statistic d40l > 



k — 



E 

4 + 



(v*-2) 
2 + q 

qg* - 1' 



(45) 



where 



V* 

i¥ 2 



and 



E" 



V* 



C2 = 



C3 



A2 
9 



1 +ciB 



(1 -c 2 B) 2 (l -c 3 B) 



(46) 



3? +2(i - 8 y 
q-g 

3q + 2(l- g y 

3<?+2(i - g y 

(q + l)Aj - (q + 4)A 2 
(q + 2)A 2 



with B,Ai,A 2 given by (J43J and d44l >. 
Section 4 in HI. 

5. Conclusions 



(47) 



For more details see 
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Appendix 

Appendix A. Derivatives of the MSE matrix with respect 
to the variance components 

Here we shall assume that G l , the inverse of G — Var(u), 
does exist, and thus we can use the MMEs as defined by |@). Al- 
though the subsequent derivation of the derivatives of the matrix 
C is general, finally we shall consider only a special case, based 
on the covariance structure of the simple linear mixed model 
( [Tol l, with the variance-covariance matrices of the following 
form: G = Var(u) - diag{cr 2 I r .}, i = 1, . . . , s, and R = Var(e) = 



Here we have presented a brief overview of the convention- 
ally used methods for making statistical inference about lin- 
ear functions of the fixed effects and/or about the fixed and 
random effects simultaneously, in conventional simple linear 
mixed model, by using the elements of the solution of the Hen- 
derson's mixed model equations. Further, we have also pre- 
sented some improvements, based on the adjusted MSE ma- 
trix of the EBLUP, as well as a generalization of the standard 
Kenward-Roger method (suggested for making statistical in- 
ference about the fixed effects) for derivation of the approxi- 
mate distribution of the Wald-type pivot statistic, suggested for 
making statistical inference about the fixed and random effects 
simultaneously. Notice that this method for derivation of the 
approximate distribution of the Wald-type pivot statistic is not 
unique. As pointed out by Alnosaier in |l|], there are several 
other alternative solutions available, however, such modifica- 
tions have not been considered here. 

The presented (explicit) expressions are valid in the simple 
LMM defined by ilOi . They are rather simple, and can be read- 
ily implemented in practically any (statistical) software envi- 
ronment. Based on the results presented in Section [Appendix ~A 
it is straightforward to get explicit expressions also for the more 
general LMM with linear parametrization of the variance-covariance 
matrices G and R, provided that the REML of variance compo- 
nents and its estimated variance-covariance matrix is available. 
The situation with nonlinear parametrization of the matrices G 
and R requires more specific approach. 



cr 2 s+l I„, soV= Var(y) = ZGZ' +R = £J=l crfZ^ + ^l n . 

Moreover, as we consider methods for statistical inference 
for estimable linear functions w = A'(b',u')' = K'b + L'u, 
i.e. such that K = X'A for some matrix A, further we shall as- 
sume, without loss of generality, that the inverse of the MME 
matrix H (the matrix on the left-hand side of the equation (|4]l) 
does exist, in particular we shall assume that the inverse of 
X'R l X does exist. Recall that 



H = (X,Z)' 'R-\X,Z) + (0,/,)'G _1 (0,/,), 



and so, 

C = H 1 or H 

Further, we shall denote 

A„ 
A, 



A J+ i 
for i = 1, 



(0,/ r )'(0,/ r ), 
(0,(0,..., I n% 

diag,{/,. 

(X,Z)'(X,Z) = H , 



0))'(0,(0,...,/ ri , 



(A.l) 

(A.2) 

(A3) 
,0)) 
(A.4) 

(A.5) 



s, where diag ( {/ f .} is (r x r)-matrix with its i-th 
diagonal block equal to I n , otherwise with zero elements. 

Further, for arbitrary matrix A we shall denote its partial 
derivatives with respect to the components of a vector parame- 
ter = (0i,...,0 i+ i)' as 

A d) _ ^ A d,j) _ d 2 A A(Ujk) _ d 3 A 
dGi' dOidef dGidBjdGk 



for i,j,k= 1 , . . . , s + 1 . 

Here we shall derive explicit expressions for derivatives of 
the matrix C, i.e. C (i) , C ( '- ;) , and C (i - j ' k \ which depend on the 
derivatives of the matrices G and R, i.e. on G w , G^, G^' k \ 
and R®, and R^. 

Recall that the derivative of A -1 , the inverse of a symmetric 
matrix A, with respect to some scalar parameter 0, is given by 



36* 



(A.7) 



and the rule for computing the derivative of a symmetric matrix 
ABA with respect to some parameter is 



8 A A A 8A OA A ifi, 
—ABA = AB— + —BA +A—A. 

80 80 80 80 



(A.8) 



Let A be an inverse of a symmetric matrix B, i.e. A = B l . 
Then, based on ( lA.7t and ( lA.8t . we define the following matrix 
operators: 

D® (A, B) = -AB (i) A, (A.9) 
D (U) (A, B) - A [B (i) AB {J) + B u) AB (i) - B iuj) ) A, (A. 1 0) 

D (.i,m (A, B) = -A (B (i) AB u) + B (j) AB (i> - B (iJ, )AB (k) A 
-AB (k) A (B (i) AB u) + B u> AB (i) - B^A 
+A (B & AB m + B (i ' k) AB ij) + B u) AB {i ' k) + B m AB (i) - 
- B (i) AB (k) AB U) _ B u) AB (k) AB (i) A + B^ k) )A. (A. 11) 



From that we directly get 

C (0 _ £>(0 (c 



(A. 12) 



= -CH (i) C, 
C (i,j) _ £)(U)( C//) 

= C (H (i) CH (j) + H u) CH (i) - H (i ' J) ) C, (A. 1 3) 
C (U* = £p-j*\c,H), (A. 14) 

for i, jf.Jt = 1, . . .,i + 1, 

Further, based on (lA.ll i. we directly get the derivatives of 



the matrix H. For i, j, k 


= l,...,s 






H® = 


(0,/,)'G 


~ m (0,I r ), 


(A.15) 




(X,Z)'/T 


Us+l \X,Z), 


(A. 16) 


H (U) = 


(0,/ f )'G 


1(U) (0,/ r ), 


(A. 17) 




(X,Z)'/T 


■ 1(S+1 > S+1 \X,Z), 


(A.18) 


H VJ,k) _ 


(0,/ r )'G 


l(iJ ' k \0,I r ), 


(A. 19) 


^-(.!+l,.!+l,S+l) 


(X,Z)'RT 


-l(.s+l,.s+l,j+l)^- 2) 


(A.20) 



where 



G -K0 = £,©( G -1 )G ) 

G -KU) _ £)(U)( G -1 ;G ) 

G -KU,*) = £)tt^)( G -l )G ) 

^-l(j+l,J+l) _ £)(.!+!, s+1) ^ 



^-l(j+l,J+l,J+l) _ £)(J+1,J+1,.!+1) ^J-l ^ (A 21) 

Notice that 



H i.um = 0) 



(A.22) 
(A.23) 



whenever one index is equal to s + 1 and some of the other 
indices is different from for s + 1 , for i,j,k=l,...,s+l. 

Appendix A.l. Derivatives of the MME matrix H in simple LMM 
In the simple LMM ( fTUI ). we get 
1 



H 



H 



(u) 



H) 



rAi, 



H) 



rAj, 



(A.24) 
(A.25) 
(A.26) 



(°f) 

for / = 1 , . . . , s + 1 . Notice that 

H {iJ) = 0, and = 0, (A.27) 

for any combination of unequal indices k = 1, . . ., s + 1. 

Appendix A.2. Derivatives of the MME matrix C in simple LMM 

By combining ( lA~T2l . ( lAJJl <E2§ , dATBl and ( TOTI 
in simple LMM ( fTUI ). we directly get 



C 



1 



c 



(i.i) 



H) 

2 



r CA,C, 



■C(A,CA,- - of a,-)c, 



(A.28) 
(A.29) 



C (U) = — I c(AiCA J + A 7 CAi)c, i* j, (A.30) 

for i, j = 1, . . . , s + 1 . 

The explicit expression for C (,,J,k \ i.e. the third partial deriva- 
tive of C for i,j,k= 1 , . . . , s + 1 , is not presented here, however, 
it can be similarly evaluated based on dA.14b . (lA.24b . ( IA.25K 
CA26I and dAT27b . 

Appendix A. 3. Derivatives of the MSE matrix M„ in simple 
LMM 

Recall that M lT ., the MSE matrix of the best linear unbiased 
predictor of w, is given by M, 7 = A'CA, where A is ((p + r)xq)- 
matrix of given coefficients. 

Let A be a solution of a system of linear equations HA = A, 
i.e. A = CA, and let A be decomposed into block-matrices such 
that A = (Ag, A'j, . . . , A' s )', where Aq is (p x ^-dimensional 
block-matrix, and A,, i = 1, are (r, x ^-dimensional 

block-matrices of A. Similarly, let {C},j denote the (i,f)-th 



bloclH of the matrix C, and let {C},. denote the z'-th row-block 
and {C}., the z'-th column-block of the matrix C. 

Then, based on the derivatives of the matrix C, we directly 
get the first partial derivatives of the MSE matrix M% with re- 
spect to the variance components <x 2 , . . . , cr}, <x 2 +1 as 



M 



1 



M' 



(s+l) 



K) 
i 

KO 2 



-A' A/A 



1 



2 A- A,-, z = 1, 



r A'A J+ iA, 



(A.31) 



where the matrices A, are defined by dA.4t and dA.51 . The sec- 
ond partial derivatives of M„ are given by: 

2 



M 



(/,/) 



(°?y 

2 

HI 



: A' (A.CA; -ct 2 A ; )a 
■ (A-{C}//A/ - <x 2 AjA,), (A.32) 



for z = 1 , . . . , s, and in for z = s + 1 we get 

2 



M 



(s+l,S+l) 



A'(A. S+1 CA J+1 -cr 2 +1 A s+1 )A,(A.33) 



Further, 



— -A' (A,CA; + A.CA,) A 

— J-j (a;{c} !7 a^ + A;.{c} 7i A«) , (a.34) 



for z ^ j, i,j — 1 , . . . , s, and 



M 



fts+i) 



= M 



(s+l,<) 



r A (A,-CA S+1 +A J+1 CA,-)A, 



- — L^(a;{C},a s+1 A 
K<i) 

+ A'A. S+1 {C},A,), (A.35) 



for i = 1 , . 



Appendix A.4. Approximation of the second component of the 
EBLUP 's MSE matrix in simple LMM 
According to ( l33l l. let us define Mg$ by 



M, 



dw 



s+l s+l 

''=1 7=1 

s+l s+l 

(=1 7=1 



3 (w - w) d(w — w) 



do-} 



(A.36) 



2 Notice that for /, j = l,...,sthe block {C)ji = {C22I1), Le. it is the (i,7)-th 
block of the matrix C22, which can be, based on 0, efficiently computed as 

C22 = cr] +l G {<r] +l I r + MG)\ where M = Z'Z - Z'X(X'X)~X'Z. 



where denote the elements of the variance-covariance matrix 
S of cr 2 . Then, by using 



w - w = A'C(X, Z)'R X (y - Xb) - A'(0, I r )'u, (A.37) 



we get 



d (w - w) 

do-} 



-A'CH (i) C(X, Z)'R-\y - Xb) 
-A'C(X,Z)'R- l R {i) R- l (y - Xb). (A.38) 



and then, by taking the covariances of the vectors with i, j = 
1, . . . , s + 1, we get, 



Cy = A'C 



H W C(X, Z)'R- l VR~ l (X, Z)CH' 



(./) 



+ H®C(X, Z)'R- 1 V/T 1 R U) R- 1 (X, Z) 
+ (X, Z)'R- l R®R- l VR- l {X, Z)CH (i) 

+ C(X, Z)'R l R (i) R- 1 VR- l R (i) R~\X, Z) 



CA, 



(A.39) 



where V = ZGZ' + R. 

Notice that in the simple LMM ([10]) we have R® = R (i) = 0, 

for i,j = l,...,s, and R (s+l) = I„. From that we get R- l R (s+l) = 
R (s+r )R -i = R -i = ij and 

Cj j = — ^A'ciAiCHvCA, 

-l{i=s+\}0- s+1 H v CAj - l(j =J+ i)cr j+1 A/C//v 

+ l (l =/= s +il (o-] +l ) 2 HyjCA, (A.40) 

for i,j = 1, . . . , s + 1, where l(,- =J+ i), ly =s +n, l{i=j =s+ i) are the 
indicator functions, and H v = {X,Z)'R x VR l {X,Z) fulfills the 
property 



CH V C = 



Cn 







G-C22 



(A.41) 



Hence, the approximation of the second component of the EBLUP's 
MSE matrix, i.e. Mg», in simple LMM is 



.s+l .s+l 
<=1 7=1 



(A.42) 



with dj, i, j = 1, . . . , s + 1, given by iAAOi . 

By recognizing that in simple LMM ( fTOt we have = 

-2C U and M^ J) = - (c, v + C ;i ), z, j = 1 5 + 1, see also 

1 22] eq. (4.6), we get the alternative expression for the approx- 
imation of the second component of the EBLUP's MSE matrix 
in simple LMM, given by 



s+l s+l 



(A.43) 



''=1 7=1 



where the matrices Mv } are given by ( lA.32b . ( lA.33b . (1A.34I I. 
and (E32). 



Appendix A.5. Bias-corrected estimator of the MSE matrix of 
EBLUP in simple LMM 

In simple LMM ( fTOb . the bias-corrected estimator of the 
MSE matrix of the empirical BLUP of w = A'(b', u')' , i.e. M A , 
is given (based on ( f39t and ( 1A.43I )). as 



Ma, + 2M 5 & 

' s+l s+l 

\i=\ j=\ 



= M-a - 



(A.44) 



and in particular, by using M a = A'CA and dA.32b . (lA.33b . 
( 1A.34I I. and ( 1A.35I ). we get 



M,, A = A' A + 



4E 5 



-A' (fr 2 s+l Ho-H CH )A 



s 4E 

- J] 7^77^2 (a;{C},^oA + A'#o{C}.A-) 



Zj Zj / „ 2 "2\ 2 



(Aj^jyAy + A^A,), (A.45) 



where A = CA, ff = A J+] = (X,Z)'(X,Z), and E, (with 
elements Ey, z',y = 1, ...,.? + 1), is the estimated variance- 
covariance matrix of the REML estimator & 2 = [& 2 , . . . , <x 2 +I ) . 
Here, A = (A ( ',, A\ , . . . , A' s ) is decomposed into block-matrices 
such that Ao is (p x g)-dimensional block-matrix, and A,, i = 
1, . . . , s, are (r,xg)-dimensional block-matrices of A. Similarly, 
{C}jj denote the (i, j')-th (r, xr/)-dimensional block of the matrix 
C, and {C},. denote the z'-th (r, x (p + r))-dimesional row-block 
and {C}., the z'-th {{p + r) x r,)-dimesional column-block of the 
matrix C. 

Appendix A.6. Generalized Kenward-Roger method for statis- 
tical inference on fixed and random effects based 
on adjusted estimator of the MSE matrix of the 
EBLUP in simple LMM 

Here we shall consider the scaled Wald-type F-statistic de- 
fined by (1401 . in particular 



kF * = — (w - w)' yil^j (yv - w) 



approx. 



(A.46) 



where M«yt is given by ( IA.45b . 

The moment based estimators of the parameters k and v are 
based on comparing the first and the second moments of the 
scaled F-statistic (1A.46E with the moments of the ^-distribution 
with q and v degrees of freedom, i.e. by solving the system of 
equations 



E{kF») =kE, = E = E(F q>v ), 
Var(KF,) = K 2 V t = V= Var{F q , v ), 



(A.47) 



where E, = E(F„) and V* = Var(F„). Based on the properties 
of the F-distribution we get 

V 

E = 



V = 



v-2 
2v 2 (v + q-2) 

?(v-2) 2 (v-4) 
2£ 2 v + q - 2 



q v-4 
provided that v > 4. By denoting 
V 



2E 2 



we get 



v = 4 + 



q + 2 
qg-l' 



(A.48) 



(A.49) 



(A.50) 



and consequently, the moment estimators of k and v are given 
as 



K = 



F,(v-2) 

a q + 1 

v = 4 + 



qg-l 



where 



V, 

Q = r. 

2£„ 2 



(A.51) 



(A.52) 



The expectation and the variance of the statistic F» defined by 
< IA.46b can be estimated by using 

E, = E(F t )=E l r-(E A (F*\fr 2 )) 
V. = Var (F.) = E & 2 (Var & (f. | <t 2 )) 

+ Var & 2 fa (F. I o- 2 )) . (A.53) 

Alnosaier in IH derived approximations for F« and V* in the 
special case, when the F-statistic ( IA.46b is restricted on fixed 
effects only. The derivation of the approximations E* and V, in 
the general case, (i.e. for the F-statistic defined by dA.461 >). is 
not presented here. However, in analogy with the derivation of 
the approximations presented in yj], we suggest F» and V», as 
the approximations of F« and V,, in the following form 

E, = 1-A 

q 

2 

V, = -(1+F), 

q 



(A.54) 



where 



B = ^(A 1+ 6A 2 ), 

s+l .s+l 

* - EE- 

'=1 /=! 
s+l s+l 

2 X 2 'V tr (M^MfM: 1 ^) . (A.55) 



E ijt r(M^M«)tr(M^MW), 



A, = 



10 



The suggested approximations depend on the unknown vari- 
ance components cr 2 = (cr 2 , . . . , cr 2 +I ) . Consequently, the sug- 
gested estimators of the parameters k and v, based on the esti- 
mated versions of ( 1A.511 I. are 



v = 4 + 



£,(v-2) 
q + 2 



qg-l 



where 



and 



e=—> 

2E„ 



E* = 1 + 



(A.56) 



(A.57) 



with 



V. = -{\+B), 



(A.58) 



At = 2Z^ tr K lM *) tr K lM *% 



s+1 j+1 



i2 = ZZ^K 1 ^* 1 ^)- (A - 59) 



In particular, by using M ff , = A'CA = A' A and ( 1A.31I I. we 
finally get 



zz 



2Ey 



i<7 



Xtrf(A'A) ' A;A,-)tr((A'A) ' A}A,) 
Xtr^A'A) 1 A;A,-V((A'A) 'A'# A) 



^tr((A'A) _I A'// A) , (A.60) 



42 - z|H( A ' A )~ Vi ) 2 ) 

+ Z Z /„ 2 -2\ 2 

xtr^(A'A) 'A;A, (a'A) 'A}A 7 j 

+z 



tr^(A'A) 'AjA^A'A) 'A'Z/oAj 

tr^(A'A)"' A'// A) 2 J, (A.61) 



x 

2jj+1,s+1 



as before, A = CA, H = A s+l = (X,Z)'(X,Z), and t, (with 
elements i, j = l,...,s + 1), is the estimated variance- 

covariance matrix of the REML estimator cr 2 - {&\, . . . , o" 2 +1 ) . 
A = Aj, . . . , A^ is decomposed into block-matrices such 
that Ao is (pxg)-dimensional block-matrix, and A,-, i = \,...,s, 
are (r, x g)-dimensional block-matrices of A. Similarly, {C},y 
denote the (z, j)-th (/*,• x r,)-dimensional block of the matrix C, 
and {C},. denote the z'-th (r, x (p + r))-dimensional row-block 
and {C}.j the j-th ((p + r) x r,)-dimensional column-block of the 
matrix C. 



Appendix B. Estimation of the variance components by solv- 
ing the MMEs 

The presented iterative procedure for estimation of the vari- 
ance components by solving the Henderson's mixed model equa- 
tions has been suggested by Searle, Casella and McCulloch in 
JH], see pp. 275-286. The MATLAB version of the algorithm 
has been implemented by Witkovsky in I68I1. 

Here we use the same notation as in [49]. In each step 
of the suggested iterative procedure, we shall denote V (,) = 
cr 2 s ( l\l r + Z'ZG {, \ G {,) = diag (of /,. ). The algorithm starts 
with the choice of the starting values for variance components 
cr 2 (°) = (cr 2(0) , . . . , of /)' and setting t = 0. In the f-th step of 
the procedure the algorithm solves the system of mixed model 
equations: 



X'X X'ZG m \( b (,) 
Z'X V {,) 



X'y 
Z'y 



(B.l) 



and B« = G w v w . 



Appendix B.l. ML estimates of the variance components 

The ML estimates of the variance components are calcu- 
lated iteratively as 



(T 



2 (f+l) 



2 C+l) 



-(f)' -W 
u. u 
i i 



i 



y (y - xh® - za w ) 

n 



(B.2) 



where a (,) is the z'-th r, -dimensional subvector of 5 W and W is 
the z'-th diagonal block of the matrix W^'\ where 

-l 



aWtf®I r + Z'Z G uy\ (B.3) 



The iterative procedure should be stopped after the f-th step 
< e, for the chosen precision limit e, and 



if Ho 2 ® -a 2 ™ 



where cr 



The final solutions of the iterative procedure are denoted by 
b, u = ffij, . . . , iV 5 ) , and cr 2 = {cr 2 , . . . , o^ +1 J ■ Similarly, we 
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denote W and use the adequate notation G, R, and C for the 
estimated versions of matrices G, R, and C. 

The log-likelihood function for ML estimation evaluated at 
the ML estimates b and <x 2 , say loglik ML , is 



loglik A 



-i«log(27r)- i log 



-\{y-Xb) \T l (y-Xb), 
= -Un log {2nfr] +i ) - log (|W|) + n) ,(B.4) 



where V = ZGZ' + ct^ +1 /„ and W = (l r + Z'ZG/<x 2 +1 ) . 

The Fisher information matrix (which is in fact the inverse 
of the asymptotic variance-covariance matrix) of the ML esti- 
mators of the variance components, say /ml(o~ 2 ), can be evalu- 
ated at the ML estimates & 2 as 



Iml(&- 2 ) = T> 



( 4 7 k,-2tr(iy, 



»„)]+tr(IV„lV j ,) \ /' 



W»'//)-EJtr(iy, ; H'j,.)V I 



ma? 

\ sW~, f. . 

Vrotv i j+1 //=! 



(B.5) 



where 6^ = 1 if (' = y, otherwise 5y = 0, and Wy is the (r,- X ry) 
block of the matrix W. 

Appendix B.2. REML estimates of the variance components 

Similarly, the REML estimates of the variance components 
are calculated iteratively as 



o-l 



C+i) 



ufuf 



; = 1 , . . . , s, 



cr 



.2 C+l) 



y' [y - Xb {,) - Zu (,) ) 



n-r x 



(B.6) 



where by r x we denote the rank of the matrix X, u. is the i- 
th r,-dimensional subvector of S w and r? is the ;-th diagonal 
block of the matrix T w , where 



(0 _ ^2 (0 



T w = cr 



s+l 



(cr} +l % + MG*f\ 



(B.7) 



where M = Z'Z - Z'X(X'X) X'Z. 

The log-likelihood function for REML estimation evaluated 
at the REML estimates cr 2 , say loglik R£ML , is 

loglik„ = -i(n-r x )log(27r)-ilog(|B't>B|) 

~y'B(B'VB) l B'y, 

= ~(n - r z )log(27ro- 2 +1 ) 

™ (- log (l?l) + («-**)), (B.8) 

where B is an nx(n-r x ) matrix, such that BB' - l n -X(X'XyX' 
and B'B = 7„_ rx . Further, f = (l r + MG/^ +1 )" . 

The Fisher information matrix of the REML estimators of 
the variance components, Ireml(<t 2 ), can be evaluated at the 
REML estimates cr 2 as 

Ireml (o -2 ) = T x 



' Mr 

I mar 



2[r(r„)]+ii-(r„r J ,- 



(row *i<^+l If 



Ji,;=l Ira/ ^+1 Ji=l 



' 



(B.9) 



where 5ij = 1 if i = j, otherwise 6/j = 0, and 7y is the (r,- x r,) 
block of the matrix 7*. 

Similarly, the final solutions of the procedure are denoted 
by h, U = ffij, . . ., ft' s ) , and cr 2 = (fr 2 , . . ., cr 2 +1 ) . Further, we 
denote T, and use the adequate notation G, R, and C for the 
estimated versions of matrices G, R, and C. 

For more details on ML and REML estimators see the Chap- 
ter 6 in Searle et al. (1992). 

Appendix B.3. MINQE's of the variance components 

For completeness, here we present procedures to calculate 
the MINQE(I) and the MINQE(U,I) estimators of the variance 
components at given (prior) values of the variance components 
2 '°? ) . Here we assume that cr 2<0) > for all 



2(0) 



1. 



2(0) 
1 ' • 

, s + 1 . For more details on minimum norm quadratic 



estimation of the variance components see e.g. [33], [42], and 

iH. 

The MINQE(I) of cr 2 , say fr 2 , at the prior value cr 2(0) is 
defined as the solution of the following system of equations 



(B.10) 



where by //(/) we denote the (s+lxs+l)-dimensional MINQE(I)- 
matrix and q - (qi,.. . ,q s +\Y denotes the vector of MINQE 
quadratic forms. The matrix is defined by its elements as 



(B.ll) 



I, 



1, where V; = Z,Z ( ', for i 

- "i i v>s+l 2i 



and y<°> = ZG^Z' - - m ' - ^ s+l - m 



V. 



s+l 



In, 



The matrix H, 



can be easily evaluated by using (IB.5t . namely 



U) 



(B.12) 



Further, the vector q of MINQE quadratic forms, defined by 
its elements as 



qi = y' (M X V {0) M X ) + Vi (M x v m) M x y 



y, 



(B.13) 



i = 1, . . . , s + 1, with M x = /„ - X(X'X)-X, could be easily 
evaluated by using 

~(0)'~(0) 

u. u 



K (0) ) 



?.s+l - 



( y _ x g(0) _ Zfi (0)y _ z g(0) _ Zfi(0 )j 



,(B.14) 



where S^ 0) is the i-th r, -dimensional subvector of S' - 1 . 

Similarly, the MINQE(U,I) of cr 2 , say cr 2 , at the prior value 
0-2(0) j s fjgfjjjgd as ); ne solution of the following system of equa- 
tions 



H(vi)& = q, 



(B.15) 
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where H^q denotes the (s+ 1X5+ l)-dimensionalMINQE(U,I) 
matrix, defined by its elements 

[H m } tJ = tr ((M x y (0) M x ) + Vi (M X V (0) M X ) + V>) , (B.16) 

i, j — 1 , . . . , s + 1 , and by using (IB. 91 ) we get 

%//)= 2 Wl (^ 2(0) ) - (B.17) 

Note that the MINQE ex 2 , defined by (lB~T0l > or by (lB~T5l 
is not given uniquely unless the MINQE matrix is of full rank. 
In fact, one version of the solution to the MINQE equations is 
<x 2 = H + q, where H + denote the Moore-Penrose ^-inverse of 
the appropriate MINQE matrix. 

The MINQE of unbiasedly estimable vector Fa 2 , where F 
is such matrix that F' = HA for some matrix A, is Fa 2 , and is 
unique. 

In particular, under given assumptions, the MINQE(U,I) 
Fa 2 , with F such that F' = H(ui)A for some matrix A, is the 
0-2(0) -locally minimum variance unbiased invariant estimator of 
Fa 2 with 

E{Fa 2 ) = Fa- 2 , 

Var{Fa 2 \a- m ) = 2FH m) F' 

= 2A'H m A. (B.18) 

On the other hand, the MINQE(I) F& 2 is a biased estimator 
of Fa 2 with 

E(F& 2 ) = FHtfI m o*, 
Var{F& 2 \a- m ) = 2FH^H (UI) H^F' . (B.19) 
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